Compare eQTLs with BXD-specific references and B6mm10 reference:
Doing this in default STAR setting (genotypesandimputed_withannotation_Local_0_10) because also available with B6 mm10 reference (abbreviated B6)
FastQTL output (http://fastqtl.sourceforge.net/):
prepareRawCounts <- function(counts, ref){
# convert gene ids to gene names
idx_uniq <- which(!duplicated(GeneConvert$V2[match(rownames(counts), GeneConvert$V1)]))
counts <- counts[idx_uniq,]
rownames(counts) <- GeneConvert$V2[match(rownames(counts), GeneConvert$V1)]
# change X into C for cortex
colnames(counts) <- gsub("X", "C", colnames(counts))
# split by tissue and condition
CortexNSD <- subset(counts, select=grep("^C[0-9]*nsd$", colnames(counts), value=TRUE))
CortexSD <- subset(counts, select=grep("^C[0-9]*sd$", colnames(counts), value=TRUE))
LiverNSD <- subset(counts, select=grep("^L[0-9]*nsd$", colnames(counts), value=TRUE))
LiverSD <- subset(counts, select=grep("^L[0-9]*sd$", colnames(counts), value=TRUE))
# clean and formate names of mouse lines
colnames(CortexNSD) <- LinesNames$V5[match(gsub("[LCnsd]", "", colnames(CortexNSD)), LinesNames$V10)]
colnames(CortexSD) <- LinesNames$V5[match(gsub("[LCnsd]", "", colnames(CortexSD)), LinesNames$V10)]
colnames(LiverNSD) <- LinesNames$V5[match(gsub("[LCnsd]", "", colnames(LiverNSD)), LinesNames$V10)]
colnames(LiverSD) <- LinesNames$V5[match(gsub("[LCnsd]", "", colnames(LiverSD)), LinesNames$V10)]
# output 4 datasets
assign(paste0("counts", ref, "ref_CortexNSD"), CortexNSD, envir=parent.frame())
assign(paste0("counts", ref, "ref_CortexSD"), CortexSD, envir=parent.frame())
assign(paste0("counts", ref, "ref_LiverNSD"), LiverNSD, envir=parent.frame())
assign(paste0("counts", ref, "ref_LiverSD"), LiverSD, envir=parent.frame())
}
# find top 3 genes for affected
top3genes <- function(tissue, condition){
t <- unlist(strsplit(tissue, split=""))[1]
dataB6 <- eval(parse(text=(paste0("dataB6", t, condition))))
dataB6full <- eval(parse(text=(paste0("dataB6full", t, condition))))
dataBXD <- eval(parse(text=(paste0("dataBXD", t, condition))))
dataBXDfull <- eval(parse(text=(paste0("dataBXDfull", t, condition))))
# identify genes in common between the B6 and BXD lists
commoneGenes <- intersect(rownames(dataB6), rownames(dataBXD))
# remove genes with 0 variants
commoneGenes <- setdiff(commoneGenes, rownames(dataB6full)[dataB6full$V2==0])
# bigger difference in slope
orderedgenes <- commoneGenes[order(abs(dataB6full[commoneGenes, "V9"]-dataBXDfull[commoneGenes, "V9"]), decreasing=TRUE)]
##dataB6full[match(orderedgenes[1:3], rownames(dataB6full)),]
##dataBXDfull[match(orderedgenes[1:3], rownames(dataBXDfull)),]
outtable <- cbind(dataB6full[orderedgenes[1:3], c(5, 8)], dataB6[orderedgenes[1:3], "adjustedpvalue"], dataBXDfull[orderedgenes[1:3], c(5, 8)], dataBXD[orderedgenes[1:3], "adjustedpvalue"])
colnames(outtable) <- c("variantB6", "slopeB6", "qvalueB6", "variantBXD", "slopeBXD", "qvalueBXD")
# output a table with info on 3 top affected genes
assign(paste0("top3genes", t, condition), outtable, envir=parent.frame())
}
# plot gene expression for a pair gene-variant
plotGeneExpressionVariant <- function(ref, gene, variant, tissue, condition){
# get first letter of tissue
t <- unlist(strsplit(tissue, split=""))[1]
# get gene expression dataset
expr <- eval(parse(text=paste0("expr", ref, "ref", t, condition)))
# plot
plot(as.numeric(expr[gene, sort(names(expr))]), pch=19, col=as.factor(genotypesGN[genotypesGN$Locus==variant, sort(names(expr))]), main=paste0("Gene: ", gene,", ref: ", ref, ", variant: ", variant, ", tissue: ", tissue, ", condition: ", condition), xaxt="n", xlab="", ylab="gene expression", las=1)
axis(1, 1:length(names(expr)), sort(names(expr)), las=2)
}
# vectorize function
plotGeneExpressionVariant <- Vectorize(plotGeneExpressionVariant, SIMPLIFY=FALSE)
# plot raw gene counts for a pair gene-variant
plotGeneCountsVariant <- function(ref, gene, variant, tissue, condition){
# get gene expression dataset
counts <- eval(parse(text=paste0("counts", ref, "ref_", tissue, condition)))
# plot
plot(as.numeric(counts[gene, sort(names(counts))]), pch=19, col=as.factor(genotypesGN[genotypesGN$Locus==variant, sort(names(counts))]), main=paste0("Gene: ", gene,", ref: ", ref, ", variant: ", variant, ", tissue: ", tissue, ", condition: ", condition), xaxt="n", xlab="", ylab="raw gene counts", las=1)
axis(1, 1:length(names(counts)), sort(names(counts)), las=2)
}
# vectorize function
plotGeneCountsVariant <- Vectorize(plotGeneCountsVariant, SIMPLIFY=FALSE)# load table to convert mouse line names
LinesNames <- read.table("F:/BXD/data/ConvertLineNames.tsv", header=FALSE, stringsAsFactors=FALSE, sep="\t")
# load table to convert between gene id and gene name
GeneConvert <- read.table("F:/BXD/references/gene_names_B6mm10.tab", header=FALSE, stringsAsFactors=FALSE, sep="\t")
# retrieve GeneNetwork (GN) genotypes
genotypesGN <- read.table("F:/BXD/data/genome/BXDGenotypes.geno", header=TRUE)
# retrieve eQTL results
for(tissue in c("Cortex", "Liver")){
t <- unlist(strsplit(tissue, split=""))[1]
for(condition in c("NSD", "SD")){
assign(paste0("dataB6", t, condition), read.table(paste0("F:/BXD/data/transcriptome/eQTL_B6mm10primaryassembly_withannotation/TMMnormalized_log2CPM_", tissue, "_", condition, "_pvalcorrected.txt"), row.names=1, stringsAsFactors=FALSE))
assign(paste0("dataB6full", t, condition), read.table(paste0("F:/BXD/data/transcriptome/eQTL_B6mm10primaryassembly_withannotation/TMMnormalized_log2CPM_", tissue, "_", condition, ".txt"), row.names=1, stringsAsFactors=FALSE))
assign(paste0("dataBXD", t, condition), read.table(paste0("F:/BXD/data/transcriptome/eQTL_BXDfromGNandimputedD2blocks_withannotation/TMMnormalized_log2CPM_", tissue, "_", condition,"_pvalcorrected.txt"), row.names=1, stringsAsFactors=FALSE))
assign(paste0("dataBXDfull", t, condition), read.table(paste0("F:/BXD/data/transcriptome/eQTL_BXDfromGNandimputedD2blocks_withannotation/TMMnormalized_log2CPM_", tissue, "_", condition, ".txt"), row.names=1, stringsAsFactors=FALSE))
}
}
# retrieve raw gene counts, formate and subset
countsB6ref <- read.table("F:/BXD/data/transcriptome/2_mapping_STAR/OnGenomeAndAnnotation/B6mm10primaryassembly_withannotation/Summary_ReadsPerGene.out.tab", header=TRUE, row.names=1)
prepareRawCounts(countsB6ref, "B6")
countsBXDref <- read.table("F:/BXD/data/transcriptome/2_mapping_STAR/OnGenomeAndAnnotation/BXDfromGNandimputedD2blocks_withannotation/Summary_ReadsPerGene.out.tab", header=TRUE, row.names=1)
prepareRawCounts(countsBXDref, "BXD")
# retrieve gene expression
exprB6refCNSD <- read.table("F:/BXD/data/transcriptome/2_mapping_STAR/OnGenomeAndAnnotation/B6mm10primaryassembly_withannotation/TMMnormalized_log2CPM_Cortex_NSD.tab")
exprBXDrefCNSD <- read.table("F:/BXD/data/transcriptome/2_mapping_STAR/OnGenomeAndAnnotation/BXDfromGNandimputedD2blocks_withannotation/TMMnormalized_log2CPM_Cortex_NSD.tab")
exprB6refCSD <- read.table("F:/BXD/data/transcriptome/2_mapping_STAR/OnGenomeAndAnnotation/B6mm10primaryassembly_withannotation/TMMnormalized_log2CPM_Cortex_SD.tab")
exprBXDrefCSD <- read.table("F:/BXD/data/transcriptome/2_mapping_STAR/OnGenomeAndAnnotation/BXDfromGNandimputedD2blocks_withannotation/TMMnormalized_log2CPM_Cortex_SD.tab")
exprB6refLNSD <- read.table("F:/BXD/data/transcriptome/2_mapping_STAR/OnGenomeAndAnnotation/B6mm10primaryassembly_withannotation/TMMnormalized_log2CPM_Liver_NSD.tab")
exprBXDrefLNSD <- read.table("F:/BXD/data/transcriptome/2_mapping_STAR/OnGenomeAndAnnotation/BXDfromGNandimputedD2blocks_withannotation/TMMnormalized_log2CPM_Liver_NSD.tab")
exprB6refLSD <- read.table("F:/BXD/data/transcriptome/2_mapping_STAR/OnGenomeAndAnnotation/B6mm10primaryassembly_withannotation/TMMnormalized_log2CPM_Liver_SD.tab")
exprBXDrefLSD <- read.table("F:/BXD/data/transcriptome/2_mapping_STAR/OnGenomeAndAnnotation/BXDfromGNandimputedD2blocks_withannotation/TMMnormalized_log2CPM_Liver_SD.tab")tissue <- "Cortex"
condition <- "NSD"
t <- unlist(strsplit(tissue, split=""))[1]
dataB6full <- eval(parse(text=paste0("dataB6full", t, condition)))
dataB6 <- eval(parse(text=paste0("dataB6", t, condition)))
dataBXDfull <- eval(parse(text=paste0("dataBXDfull", t, condition)))
dataBXD <- eval(parse(text=paste0("dataBXD", t, condition)))
# check compability of number of genes (rows)
ifelse(nrow(dataB6)==nrow(dataB6full) & nrow(dataBXD)==nrow(dataBXDfull), "Same number of genes", "Different number of genes!")## [1] "Same number of genes"
# join FastQTL statistics and q-values
eQTLB6 <- cbind(dataB6full, dataB6)
eQTLBXD <- cbind(dataBXDfull, dataBXD)
# define colnames
colnames(eQTLB6) <- c("nb_variants_tested", "MLE1_beta", "MLE2_beta", "TBD", "IDvariant", "distance", "pvalue", "slope", "pvalue_firstpermutation", "pvalue_secondpermutation", "adjustedpvalue", "marker")
colnames(eQTLBXD) <- c("nb_variants_tested", "MLE1_beta", "MLE2_beta", "TBD", "IDvariant", "distance", "pvalue", "slope", "pvalue_firstpermutation", "pvalue_secondpermutation", "adjustedpvalue", "marker")
# find common genes, not na with any reference
commongenes <- intersect(rownames(eQTLB6)[!is.na(eQTLB6[, "slope"])], rownames(eQTLBXD)[!is.na(eQTLBXD[, "slope"])])
length(commongenes)## [1] 16382
a <- eQTLB6[commongenes, "adjustedpvalue"]
b <- eQTLBXD[commongenes, "adjustedpvalue"]
diff <- abs((log10(a)-log10(b)))
plot(-log10(a), -log10(b), pch=21, bg=rgb(1,0,0, diff/max(diff)), main=paste(tissue, condition), xlab="-log10(qvalue) with B6 reference", ylab="-log10(qvalue) with BXD reference", las=1)
abline(a=0, b=1, lty=1)
abline(h=-log10(0.05), lty=2, col="darkred")
abline(v=-log10(0.05), lty=2, col="darkred")plot(-log10(a), -log10(b), pch=21, bg=rgb(1,0,0, diff/max(diff)), main=paste(tissue, condition), xlab="-log10(qvalue) with B6 reference", ylab="-log10(qvalue) with BXD reference", las=1)
text(-log10(a), -log10(b), labels=commongenes, cex=0.7, pos=3)
abline(a=0, b=1, lty=1)
abline(h=-log10(0.05), lty=2, col="darkred")
abline(v=-log10(0.05), lty=2, col="darkred")# plot slope comparison
a <- eQTLB6[commongenes, "slope"]
b <- eQTLBXD[commongenes, "slope"]
diff <- abs(a-b)
plot(a, b, pch=21, bg=rgb(1,0,0, diff/max(diff)), main=paste(tissue, condition), xlab="slope with B6 reference", ylab="slope with BXD reference", las=1)
abline(a=0, b=1, lty=1)
abline(h=0, lty=2)
abline(v=0, lty=2)plot(a, b, pch=21, bg=rgb(1,0,0, diff/max(diff)), main=paste(tissue, condition), xlab="slope with B6 reference", ylab="slope with BXD reference", las=1)
text(a, b, labels=commongenes, cex=0.7, pos=3)
abline(a=0, b=1, lty=1)
abline(h=0, lty=2)
abline(v=0, lty=2)# How many significant eQTLs are in common between the 2 references ?
sig_eQTLB6 <- commongenes[which(eQTLB6[commongenes, "adjustedpvalue"]<0.05)]
sig_eQTLBXD <- commongenes[which(eQTLBXD[commongenes, "adjustedpvalue"]<0.05)]
length(sig_eQTLB6)## [1] 4624
## [1] 4631
## [1] 4627.5
## [1] 4431
# percentage of significant eQTL in common?
length(commonsig) / mean(c(length(sig_eQTLB6), length(sig_eQTLBXD))) * 100## [1] 95.75365
# how many genes have the same marker?
changeofmarker <- commongenes[eQTLB6[commongenes, "marker"]!=eQTLBXD[commongenes, "marker"]]
length(commongenes[eQTLB6[commongenes, "marker"]==eQTLBXD[commongenes, "marker"]])## [1] 15807
## [1] 16382
# percentage of genes with the same marker?
length(commongenes[eQTLB6[commongenes, "marker"]==eQTLBXD[commongenes, "marker"]]) / length(commongenes) * 100## [1] 96.49005
## [1] 4334
## [1] 93.65748
##intersect(commonsig, changeofmarker)
##intersect(sig_eQTLB6, changeofmarker)
##intersect(sig_eQTLBXD, changeofmarker)
# Is this marginal differences (around FDR 0.05)
uniqsigB6 <- setdiff(sig_eQTLB6, sig_eQTLBXD)
length(uniqsigB6)## [1] 193
## [1] 200
hist(-log10(eQTLB6[uniqsigB6, "adjustedpvalue"]), breaks=50, las=1)
abline(v=-log10(0.05), lty=2, col="darkred")hist(-log10(eQTLBXD[uniqsigBXD, "adjustedpvalue"]), breaks=50, las=1)
abline(v=-log10(0.05), lty=2, col="darkred")# eQTL with same sign of slope
changeofsign <- commongenes[sign(eQTLB6[commongenes, "slope"])!=sign(eQTLBXD[commongenes, "slope"])]
##intersect(commonsig, changeofsign)
##intersect(commonsig, changeofmarker)
##intersect(changeofsign, changeofmarker)
##setdiff(setdiff(commonsig, changeofsign), changeofmarker)
# significant eQTLs maintained
##setdiff(setdiff(commonsig, changeofsign), changeofmarker)
length(setdiff(setdiff(commonsig, changeofsign), changeofmarker))## [1] 4323
## [1] 4431
## [1] 97.56263
# eQTLs maintained
##setdiff(setdiff(commongenes, changeofsign), changeofmarker)
length(setdiff(setdiff(commongenes, changeofsign), changeofmarker))## [1] 15745
## [1] 16382
## [1] 96.11159
a <- eQTLB6[commongenes, "adjustedpvalue"]
b <- eQTLB6[commongenes, "slope"]
plot(-log10(a), b, pch=21, bg=rgb(0,0,0, 0.1), main=paste(tissue, condition), xlab="-log10(qvalue) with B6 reference", ylab="slope with B6 reference", las=1)
abline(h=0, lty=1)
abline(v=-log10(0.05), lty=2, col="darkred")plot(-log10(a), b, pch=21, bg=rgb(0,0,0, 0.1), main=paste(tissue, condition), xlab="-log10(qvalue) with B6 reference", ylab="slope with B6 reference", las=1)
abline(h=0, lty=1)
abline(v=-log10(0.05), lty=2, col="darkred")
text(-log10(a), b, labels=commongenes, cex=0.7, pos=3)a <- eQTLBXD[commongenes, "adjustedpvalue"]
b <- eQTLBXD[commongenes, "slope"]
plot(-log10(a), b, pch=21, bg=rgb(0,0,0, 0.1), main=paste(tissue, condition), xlab="-log10(qvalue) with BXD reference", ylab="slope with BXD reference", las=1)
abline(h=0, lty=1)
abline(v=-log10(0.05), lty=2, col="darkred")plot(-log10(a), b, pch=21, bg=rgb(0,0,0, 0.1), main=paste(tissue, condition), xlab="-log10(qvalue) with BXD reference", ylab="slope with BXD reference", las=1)
abline(h=0, lty=1)
abline(v=-log10(0.05), lty=2, col="darkred")
text(-log10(a), b, labels=commongenes, cex=0.7, pos=3)# plot diagnostic for p-values
B6 <- eQTLB6[, "pvalue"]
BXD <- eQTLBXD[, "pvalue"]
hist(B6, breaks=50, col=rgb(1,0,0,0.5), xlab="p-values",
ylab="frequency", main=paste("Diagnostic p-values plot", tissue, condition), las=1)
hist(BXD, breaks=50, col=rgb(0,0,1,0.5), add=TRUE)
# add legend
legend("topright", legend=c("B6","BXD"), col=c(rgb(1,0,0,0.5), rgb(0,0,1,0.5)), pt.cex=2, pch=15)# plot diagnostic for q-values
B6 <- eQTLB6[, "adjustedpvalue"]
BXD <- eQTLBXD[, "adjustedpvalue"]
hist(B6, breaks=50, col=rgb(1,0,0,0.5), xlab="adjusted p-values",
ylab="frequency", main=paste("Diagnostic adjusted p-values plot", tissue, condition), las=1)
hist(BXD, breaks=50, col=rgb(0,0,1,0.5), add=TRUE)
# add legend
legend("topright", legend=c("B6","BXD"), col=c(rgb(1,0,0,0.5), rgb(0,0,1,0.5)), pt.cex=2, pch=15)# percentage of significant eQTLs over expressed genes
length(which(eQTLB6[, "adjustedpvalue"]<0.05)) / length(eQTLB6[, "adjustedpvalue"]) * 100## [1] 28.07294
## [1] 28.22664
qB6 <- qvalue::qvalue(eQTLB6$pvalue_secondpermutation)
qBXD <- qvalue::qvalue(eQTLBXD$pvalue_secondpermutation)
summary(qB6)##
## Call:
## qvalue::qvalue(p = eQTLB6$pvalue_secondpermutation)
##
## pi0: 0.5596393
##
## Cumulative number of significant calls:
##
## <1e-04 <0.001 <0.01 <0.025 <0.05 <0.1 <1
## p-value 2376 2988 3950 4629 5298 6321 16452
## q-value 2067 2614 3466 4000 4634 5539 16452
## local FDR 1668 2041 2613 2935 3231 3635 16440
##
## Call:
## qvalue::qvalue(p = eQTLBXD$pvalue_secondpermutation)
##
## pi0: 0.5580651
##
## Cumulative number of significant calls:
##
## <1e-04 <0.001 <0.01 <0.025 <0.05 <0.1 <1
## p-value 2353 2957 3950 4629 5276 6282 16395
## q-value 2061 2594 3463 3999 4643 5501 16395
## local FDR 1635 2031 2599 2922 3250 3640 16379
tissue <- "Cortex"
condition <- "SD"
t <- unlist(strsplit(tissue, split=""))[1]
dataB6full <- eval(parse(text=paste0("dataB6full", t, condition)))
dataB6 <- eval(parse(text=paste0("dataB6", t, condition)))
dataBXDfull <- eval(parse(text=paste0("dataBXDfull", t, condition)))
dataBXD <- eval(parse(text=paste0("dataBXD", t, condition)))
# check compability of number of genes (rows)
ifelse(nrow(dataB6)==nrow(dataB6full) & nrow(dataBXD)==nrow(dataBXDfull), "Same number of genes", "Different number of genes!")## [1] "Same number of genes"
# join FastQTL statistics and q-values
eQTLB6 <- cbind(dataB6full, dataB6)
eQTLBXD <- cbind(dataBXDfull, dataBXD)
# define colnames
colnames(eQTLB6) <- c("nb_variants_tested", "MLE1_beta", "MLE2_beta", "TBD", "IDvariant", "distance", "pvalue", "slope", "pvalue_firstpermutation", "pvalue_secondpermutation", "adjustedpvalue", "marker")
colnames(eQTLBXD) <- c("nb_variants_tested", "MLE1_beta", "MLE2_beta", "TBD", "IDvariant", "distance", "pvalue", "slope", "pvalue_firstpermutation", "pvalue_secondpermutation", "adjustedpvalue", "marker")
# find common genes, not na with any reference
commongenes <- intersect(rownames(eQTLB6)[!is.na(eQTLB6[, "slope"])], rownames(eQTLBXD)[!is.na(eQTLBXD[, "slope"])])
length(commongenes)## [1] 16382
a <- eQTLB6[commongenes, "adjustedpvalue"]
b <- eQTLBXD[commongenes, "adjustedpvalue"]
diff <- abs((log10(a)-log10(b)))
plot(-log10(a), -log10(b), pch=21, bg=rgb(1,0,0, diff/max(diff)), main=paste(tissue, condition), xlab="-log10(qvalue) with B6 reference", ylab="-log10(qvalue) with BXD reference", las=1)
abline(a=0, b=1, lty=1)
abline(h=-log10(0.05), lty=2, col="darkred")
abline(v=-log10(0.05), lty=2, col="darkred")plot(-log10(a), -log10(b), pch=21, bg=rgb(1,0,0, diff/max(diff)), main=paste(tissue, condition), xlab="-log10(qvalue) with B6 reference", ylab="-log10(qvalue) with BXD reference", las=1)
text(-log10(a), -log10(b), labels=commongenes, cex=0.7, pos=3)
abline(a=0, b=1, lty=1)
abline(h=-log10(0.05), lty=2, col="darkred")
abline(v=-log10(0.05), lty=2, col="darkred")# plot slope comparison
a <- eQTLB6[commongenes, "slope"]
b <- eQTLBXD[commongenes, "slope"]
diff <- abs(a-b)
plot(a, b, pch=21, bg=rgb(1,0,0, diff/max(diff)), main=paste(tissue, condition), xlab="slope with B6 reference", ylab="slope with BXD reference", las=1)
abline(a=0, b=1, lty=1)
abline(h=0, lty=2)
abline(v=0, lty=2)plot(a, b, pch=21, bg=rgb(1,0,0, diff/max(diff)), main=paste(tissue, condition), xlab="slope with B6 reference", ylab="slope with BXD reference", las=1)
text(a, b, labels=commongenes, cex=0.7, pos=3)
abline(a=0, b=1, lty=1)
abline(h=0, lty=2)
abline(v=0, lty=2)# How many significant eQTLs are in common between the 2 references ?
sig_eQTLB6 <- commongenes[which(eQTLB6[commongenes, "adjustedpvalue"]<0.05)]
sig_eQTLBXD <- commongenes[which(eQTLBXD[commongenes, "adjustedpvalue"]<0.05)]
length(sig_eQTLB6)## [1] 5018
## [1] 5020
## [1] 5019
## [1] 4834
# percentage of significant eQTL in common?
length(commonsig) / mean(c(length(sig_eQTLB6), length(sig_eQTLBXD))) * 100## [1] 96.31401
# how many genes have the same marker?
changeofmarker <- commongenes[eQTLB6[commongenes, "marker"]!=eQTLBXD[commongenes, "marker"]]
length(commongenes[eQTLB6[commongenes, "marker"]==eQTLBXD[commongenes, "marker"]])## [1] 15729
## [1] 16382
# percentage of genes with the same marker?
length(commongenes[eQTLB6[commongenes, "marker"]==eQTLBXD[commongenes, "marker"]]) / length(commongenes) * 100## [1] 96.01392
## [1] 4727
## [1] 94.18211
##intersect(commonsig, changeofmarker)
##intersect(sig_eQTLB6, changeofmarker)
##intersect(sig_eQTLBXD, changeofmarker)
# Is this marginal differences (around FDR 0.05)
uniqsigB6 <- setdiff(sig_eQTLB6, sig_eQTLBXD)
length(uniqsigB6)## [1] 184
## [1] 186
hist(-log10(eQTLB6[uniqsigB6, "adjustedpvalue"]), breaks=50, las=1)
abline(v=-log10(0.05), lty=2, col="darkred")hist(-log10(eQTLBXD[uniqsigBXD, "adjustedpvalue"]), breaks=50, las=1)
abline(v=-log10(0.05), lty=2, col="darkred")# eQTL with same sign of slope
changeofsign <- commongenes[sign(eQTLB6[commongenes, "slope"])!=sign(eQTLBXD[commongenes, "slope"])]
##intersect(commonsig, changeofsign)
##intersect(commonsig, changeofmarker)
##intersect(changeofsign, changeofmarker)
##setdiff(setdiff(commonsig, changeofsign), changeofmarker)
# significant eQTLs maintained
##setdiff(setdiff(commonsig, changeofsign), changeofmarker)
length(setdiff(setdiff(commonsig, changeofsign), changeofmarker))## [1] 4714
## [1] 4834
## [1] 97.51758
# eQTLs maintained
##setdiff(setdiff(commongenes, changeofsign), changeofmarker)
length(setdiff(setdiff(commongenes, changeofsign), changeofmarker))## [1] 15667
## [1] 16382
## [1] 95.63545
a <- eQTLB6[commongenes, "adjustedpvalue"]
b <- eQTLB6[commongenes, "slope"]
plot(-log10(a), b, pch=21, bg=rgb(0,0,0, 0.1), main=paste(tissue, condition), xlab="-log10(qvalue) with B6 reference", ylab="slope with B6 reference", las=1)
abline(h=0, lty=1)
abline(v=-log10(0.05), lty=2, col="darkred")plot(-log10(a), b, pch=21, bg=rgb(0,0,0, 0.1), main=paste(tissue, condition), xlab="-log10(qvalue) with B6 reference", ylab="slope with B6 reference", las=1)
abline(h=0, lty=1)
abline(v=-log10(0.05), lty=2, col="darkred")
text(-log10(a), b, labels=commongenes, cex=0.7, pos=3)a <- eQTLBXD[commongenes, "adjustedpvalue"]
b <- eQTLBXD[commongenes, "slope"]
plot(-log10(a), b, pch=21, bg=rgb(0,0,0, 0.1), main=paste(tissue, condition), xlab="-log10(qvalue) with BXD reference", ylab="slope with BXD reference", las=1)
abline(h=0, lty=1)
abline(v=-log10(0.05), lty=2, col="darkred")plot(-log10(a), b, pch=21, bg=rgb(0,0,0, 0.1), main=paste(tissue, condition), xlab="-log10(qvalue) with BXD reference", ylab="slope with BXD reference", las=1)
abline(h=0, lty=1)
abline(v=-log10(0.05), lty=2, col="darkred")
text(-log10(a), b, labels=commongenes, cex=0.7, pos=3)# plot diagnostic for p-values
B6 <- eQTLB6[, "pvalue"]
BXD <- eQTLBXD[, "pvalue"]
hist(B6, breaks=50, col=rgb(1,0,0,0.5), xlab="p-values",
ylab="frequency", main=paste("Diagnostic p-values plot", tissue, condition), las=1)
hist(BXD, breaks=50, col=rgb(0,0,1,0.5), add=TRUE)
# add legend
legend("topright", legend=c("B6","BXD"), col=c(rgb(1,0,0,0.5), rgb(0,0,1,0.5)), pt.cex=2, pch=15)# plot diagnostic for q-values
B6 <- eQTLB6[, "adjustedpvalue"]
BXD <- eQTLBXD[, "adjustedpvalue"]
hist(B6, breaks=50, col=rgb(1,0,0,0.5), xlab="adjusted p-values",
ylab="frequency", main=paste("Diagnostic adjusted p-values plot", tissue, condition), las=1)
hist(BXD, breaks=50, col=rgb(0,0,1,0.5), add=TRUE)
# add legend
legend("topright", legend=c("B6","BXD"), col=c(rgb(1,0,0,0.5), rgb(0,0,1,0.5)), pt.cex=2, pch=15)# percentage of significant eQTLs over expressed genes
length(which(eQTLB6[, "adjustedpvalue"]<0.05)) / length(eQTLB6[, "adjustedpvalue"]) * 100## [1] 30.4598
## [1] 30.58545
qB6 <- qvalue::qvalue(eQTLB6$pvalue_secondpermutation)
qBXD <- qvalue::qvalue(eQTLBXD$pvalue_secondpermutation)
summary(qB6)##
## Call:
## qvalue::qvalue(p = eQTLB6$pvalue_secondpermutation)
##
## pi0: 0.5149926
##
## Cumulative number of significant calls:
##
## <1e-04 <0.001 <0.01 <0.025 <0.05 <0.1 <1
## p-value 2569 3195 4224 4882 5627 6577 16451
## q-value 2278 2862 3775 4398 5028 6054 16452
## local FDR 1812 2243 2853 3204 3518 3985 16443
##
## Call:
## qvalue::qvalue(p = eQTLBXD$pvalue_secondpermutation)
##
## pi0: 0.5160343
##
## Cumulative number of significant calls:
##
## <1e-04 <0.001 <0.01 <0.025 <0.05 <0.1 <1
## p-value 2540 3196 4218 4887 5545 6555 16395
## q-value 2265 2852 3766 4411 5031 6012 16395
## local FDR 1798 2233 2851 3210 3526 3973 16388
tissue <- "Liver"
condition <- "NSD"
t <- unlist(strsplit(tissue, split=""))[1]
dataB6full <- eval(parse(text=paste0("dataB6full", t, condition)))
dataB6 <- eval(parse(text=paste0("dataB6", t, condition)))
dataBXDfull <- eval(parse(text=paste0("dataBXDfull", t, condition)))
dataBXD <- eval(parse(text=paste0("dataBXD", t, condition)))
# check compability of number of genes (rows)
ifelse(nrow(dataB6)==nrow(dataB6full) & nrow(dataBXD)==nrow(dataBXDfull), "Same number of genes", "Different number of genes!")## [1] "Same number of genes"
# join FastQTL statistics and q-values
eQTLB6 <- cbind(dataB6full, dataB6)
eQTLBXD <- cbind(dataBXDfull, dataBXD)
# define colnames
colnames(eQTLB6) <- c("nb_variants_tested", "MLE1_beta", "MLE2_beta", "TBD", "IDvariant", "distance", "pvalue", "slope", "pvalue_firstpermutation", "pvalue_secondpermutation", "adjustedpvalue", "marker")
colnames(eQTLBXD) <- c("nb_variants_tested", "MLE1_beta", "MLE2_beta", "TBD", "IDvariant", "distance", "pvalue", "slope", "pvalue_firstpermutation", "pvalue_secondpermutation", "adjustedpvalue", "marker")
# find common genes, not na with any reference
commongenes <- intersect(rownames(eQTLB6)[!is.na(eQTLB6[, "slope"])], rownames(eQTLBXD)[!is.na(eQTLBXD[, "slope"])])
length(commongenes)## [1] 12752
a <- eQTLB6[commongenes, "adjustedpvalue"]
b <- eQTLBXD[commongenes, "adjustedpvalue"]
diff <- abs((log10(a)-log10(b)))
plot(-log10(a), -log10(b), pch=21, bg=rgb(1,0,0, diff/max(diff)), main=paste(tissue, condition), xlab="-log10(qvalue) with B6 reference", ylab="-log10(qvalue) with BXD reference", las=1)
abline(a=0, b=1, lty=1)
abline(h=-log10(0.05), lty=2, col="darkred")
abline(v=-log10(0.05), lty=2, col="darkred")plot(-log10(a), -log10(b), pch=21, bg=rgb(1,0,0, diff/max(diff)), main=paste(tissue, condition), xlab="-log10(qvalue) with B6 reference", ylab="-log10(qvalue) with BXD reference", las=1)
text(-log10(a), -log10(b), labels=commongenes, cex=0.7, pos=3)
abline(a=0, b=1, lty=1)
abline(h=-log10(0.05), lty=2, col="darkred")
abline(v=-log10(0.05), lty=2, col="darkred")# plot slope comparison
a <- eQTLB6[commongenes, "slope"]
b <- eQTLBXD[commongenes, "slope"]
diff <- abs(a-b)
plot(a, b, pch=21, bg=rgb(1,0,0, diff/max(diff)), main=paste(tissue, condition), xlab="slope with B6 reference", ylab="slope with BXD reference", las=1)
abline(a=0, b=1, lty=1)
abline(h=0, lty=2)
abline(v=0, lty=2)plot(a, b, pch=21, bg=rgb(1,0,0, diff/max(diff)), main=paste(tissue, condition), xlab="slope with B6 reference", ylab="slope with BXD reference", las=1)
text(a, b, labels=commongenes, cex=0.7, pos=3)
abline(a=0, b=1, lty=1)
abline(h=0, lty=2)
abline(v=0, lty=2)# How many significant eQTLs are in common between the 2 references ?
sig_eQTLB6 <- commongenes[which(eQTLB6[commongenes, "adjustedpvalue"]<0.05)]
sig_eQTLBXD <- commongenes[which(eQTLBXD[commongenes, "adjustedpvalue"]<0.05)]
length(sig_eQTLB6)## [1] 3192
## [1] 3202
## [1] 3197
## [1] 3077
# percentage of significant eQTL in common?
length(commonsig) / mean(c(length(sig_eQTLB6), length(sig_eQTLBXD))) * 100## [1] 96.24648
# how many genes have the same marker?
changeofmarker <- commongenes[eQTLB6[commongenes, "marker"]!=eQTLBXD[commongenes, "marker"]]
length(commongenes[eQTLB6[commongenes, "marker"]==eQTLBXD[commongenes, "marker"]])## [1] 12467
## [1] 12752
# percentage of genes with the same marker?
length(commongenes[eQTLB6[commongenes, "marker"]==eQTLBXD[commongenes, "marker"]]) / length(commongenes) * 100## [1] 97.76506
## [1] 3032
## [1] 94.83891
##intersect(commonsig, changeofmarker)
##intersect(sig_eQTLB6, changeofmarker)
##intersect(sig_eQTLBXD, changeofmarker)
# Is this marginal differences (around FDR 0.05)
uniqsigB6 <- setdiff(sig_eQTLB6, sig_eQTLBXD)
length(uniqsigB6)## [1] 115
## [1] 125
hist(-log10(eQTLB6[uniqsigB6, "adjustedpvalue"]), breaks=50, las=1)
abline(v=-log10(0.05), lty=2, col="darkred")hist(-log10(eQTLBXD[uniqsigBXD, "adjustedpvalue"]), breaks=50, las=1)
abline(v=-log10(0.05), lty=2, col="darkred")# eQTL with same sign of slope
changeofsign <- commongenes[sign(eQTLB6[commongenes, "slope"])!=sign(eQTLBXD[commongenes, "slope"])]
##intersect(commonsig, changeofsign)
##intersect(commonsig, changeofmarker)
##intersect(changeofsign, changeofmarker)
##setdiff(setdiff(commonsig, changeofsign), changeofmarker)
# significant eQTLs maintained
##setdiff(setdiff(commonsig, changeofsign), changeofmarker)
length(setdiff(setdiff(commonsig, changeofsign), changeofmarker))## [1] 3026
## [1] 3077
## [1] 98.34254
# eQTLs maintained
##setdiff(setdiff(commongenes, changeofsign), changeofmarker)
length(setdiff(setdiff(commongenes, changeofsign), changeofmarker))## [1] 12435
## [1] 12752
## [1] 97.51412
a <- eQTLB6[commongenes, "adjustedpvalue"]
b <- eQTLB6[commongenes, "slope"]
plot(-log10(a), b, pch=21, bg=rgb(0,0,0, 0.1), main=paste(tissue, condition), xlab="-log10(qvalue) with B6 reference", ylab="slope with B6 reference", las=1)
abline(h=0, lty=1)
abline(v=-log10(0.05), lty=2, col="darkred")plot(-log10(a), b, pch=21, bg=rgb(0,0,0, 0.1), main=paste(tissue, condition), xlab="-log10(qvalue) with B6 reference", ylab="slope with B6 reference", las=1)
abline(h=0, lty=1)
abline(v=-log10(0.05), lty=2, col="darkred")
text(-log10(a), b, labels=commongenes, cex=0.7, pos=3)a <- eQTLBXD[commongenes, "adjustedpvalue"]
b <- eQTLBXD[commongenes, "slope"]
plot(-log10(a), b, pch=21, bg=rgb(0,0,0, 0.1), main=paste(tissue, condition), xlab="-log10(qvalue) with BXD reference", ylab="slope with BXD reference", las=1)
abline(h=0, lty=1)
abline(v=-log10(0.05), lty=2, col="darkred")plot(-log10(a), b, pch=21, bg=rgb(0,0,0, 0.1), main=paste(tissue, condition), xlab="-log10(qvalue) with BXD reference", ylab="slope with BXD reference", las=1)
abline(h=0, lty=1)
abline(v=-log10(0.05), lty=2, col="darkred")
text(-log10(a), b, labels=commongenes, cex=0.7, pos=3)# plot diagnostic for p-values
B6 <- eQTLB6[, "pvalue"]
BXD <- eQTLBXD[, "pvalue"]
hist(B6, breaks=50, col=rgb(1,0,0,0.5), xlab="p-values",
ylab="frequency", main=paste("Diagnostic p-values plot", tissue, condition), las=1)
hist(BXD, breaks=50, col=rgb(0,0,1,0.5), add=TRUE)
# add legend
legend("topright", legend=c("B6","BXD"), col=c(rgb(1,0,0,0.5), rgb(0,0,1,0.5)), pt.cex=2, pch=15)# plot diagnostic for q-values
B6 <- eQTLB6[, "adjustedpvalue"]
BXD <- eQTLBXD[, "adjustedpvalue"]
hist(B6, breaks=50, col=rgb(1,0,0,0.5), xlab="adjusted p-values",
ylab="frequency", main=paste("Diagnostic adjusted p-values plot", tissue, condition), las=1)
hist(BXD, breaks=50, col=rgb(0,0,1,0.5), add=TRUE)
# add legend
legend("topright", legend=c("B6","BXD"), col=c(rgb(1,0,0,0.5), rgb(0,0,1,0.5)), pt.cex=2, pch=15)# percentage of significant eQTLs over expressed genes
length(which(eQTLB6[, "adjustedpvalue"]<0.05)) / length(eQTLB6[, "adjustedpvalue"]) * 100## [1] 24.97083
## [1] 25.09767
qB6 <- qvalue::qvalue(eQTLB6$pvalue_secondpermutation)
qBXD <- qvalue::qvalue(eQTLBXD$pvalue_secondpermutation)
summary(qB6)##
## Call:
## qvalue::qvalue(p = eQTLB6$pvalue_secondpermutation)
##
## pi0: 0.5610873
##
## Cumulative number of significant calls:
##
## <1e-04 <0.001 <0.01 <0.025 <0.05 <0.1 <1
## p-value 1601 2031 2803 3289 3842 4647 12819
## q-value 1368 1743 2347 2783 3210 3905 12819
## local FDR 1102 1356 1751 1983 2207 2531 12795
##
## Call:
## qvalue::qvalue(p = eQTLBXD$pvalue_secondpermutation)
##
## pi0: 0.5532508
##
## Cumulative number of significant calls:
##
## <1e-04 <0.001 <0.01 <0.025 <0.05 <0.1 <1
## p-value 1589 2019 2782 3283 3833 4642 12764
## q-value 1351 1735 2345 2775 3212 3941 12764
## local FDR 1074 1339 1744 1977 2202 2538 12739
tissue <- "Liver"
condition <- "SD"
t <- unlist(strsplit(tissue, split=""))[1]
dataB6full <- eval(parse(text=paste0("dataB6full", t, condition)))
dataB6 <- eval(parse(text=paste0("dataB6", t, condition)))
dataBXDfull <- eval(parse(text=paste0("dataBXDfull", t, condition)))
dataBXD <- eval(parse(text=paste0("dataBXD", t, condition)))
# check compability of number of genes (rows)
ifelse(nrow(dataB6)==nrow(dataB6full) & nrow(dataBXD)==nrow(dataBXDfull), "Same number of genes", "Different number of genes!")## [1] "Same number of genes"
# join FastQTL statistics and q-values
eQTLB6 <- cbind(dataB6full, dataB6)
eQTLBXD <- cbind(dataBXDfull, dataBXD)
# define colnames
colnames(eQTLB6) <- c("nb_variants_tested", "MLE1_beta", "MLE2_beta", "TBD", "IDvariant", "distance", "pvalue", "slope", "pvalue_firstpermutation", "pvalue_secondpermutation", "adjustedpvalue", "marker")
colnames(eQTLBXD) <- c("nb_variants_tested", "MLE1_beta", "MLE2_beta", "TBD", "IDvariant", "distance", "pvalue", "slope", "pvalue_firstpermutation", "pvalue_secondpermutation", "adjustedpvalue", "marker")
# find common genes, not na with any reference
commongenes <- intersect(rownames(eQTLB6)[!is.na(eQTLB6[, "slope"])], rownames(eQTLBXD)[!is.na(eQTLBXD[, "slope"])])
length(commongenes)## [1] 12752
a <- eQTLB6[commongenes, "adjustedpvalue"]
b <- eQTLBXD[commongenes, "adjustedpvalue"]
diff <- abs((log10(a)-log10(b)))
plot(-log10(a), -log10(b), pch=21, bg=rgb(1,0,0, diff/max(diff)), main=paste(tissue, condition), xlab="-log10(qvalue) with B6 reference", ylab="-log10(qvalue) with BXD reference", las=1)
abline(a=0, b=1, lty=1)
abline(h=-log10(0.05), lty=2, col="darkred")
abline(v=-log10(0.05), lty=2, col="darkred")plot(-log10(a), -log10(b), pch=21, bg=rgb(1,0,0, diff/max(diff)), main=paste(tissue, condition), xlab="-log10(qvalue) with B6 reference", ylab="-log10(qvalue) with BXD reference", las=1)
text(-log10(a), -log10(b), labels=commongenes, cex=0.7, pos=3)
abline(a=0, b=1, lty=1)
abline(h=-log10(0.05), lty=2, col="darkred")
abline(v=-log10(0.05), lty=2, col="darkred")# plot slope comparison
a <- eQTLB6[commongenes, "slope"]
b <- eQTLBXD[commongenes, "slope"]
diff <- abs(a-b)
plot(a, b, pch=21, bg=rgb(1,0,0, diff/max(diff)), main=paste(tissue, condition), xlab="slope with B6 reference", ylab="slope with BXD reference", las=1)
abline(a=0, b=1, lty=1)
abline(h=0, lty=2)
abline(v=0, lty=2)plot(a, b, pch=21, bg=rgb(1,0,0, diff/max(diff)), main=paste(tissue, condition), xlab="slope with B6 reference", ylab="slope with BXD reference", las=1)
text(a, b, labels=commongenes, cex=0.7, pos=3)
abline(a=0, b=1, lty=1)
abline(h=0, lty=2)
abline(v=0, lty=2)# How many significant eQTLs are in common between the 2 references ?
sig_eQTLB6 <- commongenes[which(eQTLB6[commongenes, "adjustedpvalue"]<0.05)]
sig_eQTLBXD <- commongenes[which(eQTLBXD[commongenes, "adjustedpvalue"]<0.05)]
length(sig_eQTLB6)## [1] 2817
## [1] 2838
## [1] 2827.5
## [1] 2712
# percentage of significant eQTL in common?
length(commonsig) / mean(c(length(sig_eQTLB6), length(sig_eQTLBXD))) * 100## [1] 95.91512
# how many genes have the same marker?
changeofmarker <- commongenes[eQTLB6[commongenes, "marker"]!=eQTLBXD[commongenes, "marker"]]
length(commongenes[eQTLB6[commongenes, "marker"]==eQTLBXD[commongenes, "marker"]])## [1] 12470
## [1] 12752
# percentage of genes with the same marker?
length(commongenes[eQTLB6[commongenes, "marker"]==eQTLBXD[commongenes, "marker"]]) / length(commongenes) * 100## [1] 97.78858
## [1] 2681
## [1] 94.81874
##intersect(commonsig, changeofmarker)
##intersect(sig_eQTLB6, changeofmarker)
##intersect(sig_eQTLBXD, changeofmarker)
# Is this marginal differences (around FDR 0.05)
uniqsigB6 <- setdiff(sig_eQTLB6, sig_eQTLBXD)
length(uniqsigB6)## [1] 105
## [1] 126
hist(-log10(eQTLB6[uniqsigB6, "adjustedpvalue"]), breaks=50, las=1)
abline(v=-log10(0.05), lty=2, col="darkred")hist(-log10(eQTLBXD[uniqsigBXD, "adjustedpvalue"]), breaks=50, las=1)
abline(v=-log10(0.05), lty=2, col="darkred")# eQTL with same sign of slope
changeofsign <- commongenes[sign(eQTLB6[commongenes, "slope"])!=sign(eQTLBXD[commongenes, "slope"])]
##intersect(commonsig, changeofsign)
##intersect(commonsig, changeofmarker)
##intersect(changeofsign, changeofmarker)
##setdiff(setdiff(commonsig, changeofsign), changeofmarker)
# significant eQTLs maintained
##setdiff(setdiff(commonsig, changeofsign), changeofmarker)
length(setdiff(setdiff(commonsig, changeofsign), changeofmarker))## [1] 2678
## [1] 2712
## [1] 98.74631
# eQTLs maintained
##setdiff(setdiff(commongenes, changeofsign), changeofmarker)
length(setdiff(setdiff(commongenes, changeofsign), changeofmarker))## [1] 12444
## [1] 12752
## [1] 97.58469
a <- eQTLB6[commongenes, "adjustedpvalue"]
b <- eQTLB6[commongenes, "slope"]
plot(-log10(a), b, pch=21, bg=rgb(0,0,0, 0.1), main=paste(tissue, condition), xlab="-log10(qvalue) with B6 reference", ylab="slope with B6 reference", las=1)
abline(h=0, lty=1)
abline(v=-log10(0.05), lty=2, col="darkred")plot(-log10(a), b, pch=21, bg=rgb(0,0,0, 0.1), main=paste(tissue, condition), xlab="-log10(qvalue) with B6 reference", ylab="slope with B6 reference", las=1)
abline(h=0, lty=1)
abline(v=-log10(0.05), lty=2, col="darkred")
text(-log10(a), b, labels=commongenes, cex=0.7, pos=3)a <- eQTLBXD[commongenes, "adjustedpvalue"]
b <- eQTLBXD[commongenes, "slope"]
plot(-log10(a), b, pch=21, bg=rgb(0,0,0, 0.1), main=paste(tissue, condition), xlab="-log10(qvalue) with BXD reference", ylab="slope with BXD reference", las=1)
abline(h=0, lty=1)
abline(v=-log10(0.05), lty=2, col="darkred")plot(-log10(a), b, pch=21, bg=rgb(0,0,0, 0.1), main=paste(tissue, condition), xlab="-log10(qvalue) with BXD reference", ylab="slope with BXD reference", las=1)
abline(h=0, lty=1)
abline(v=-log10(0.05), lty=2, col="darkred")
text(-log10(a), b, labels=commongenes, cex=0.7, pos=3)# plot diagnostic for p-values
B6 <- eQTLB6[, "pvalue"]
BXD <- eQTLBXD[, "pvalue"]
hist(B6, breaks=50, col=rgb(1,0,0,0.5), xlab="p-values",
ylab="frequency", main=paste("Diagnostic p-values plot", tissue, condition), las=1)
hist(BXD, breaks=50, col=rgb(0,0,1,0.5), add=TRUE)
# add legend
legend("topright", legend=c("B6","BXD"), col=c(rgb(1,0,0,0.5), rgb(0,0,1,0.5)), pt.cex=2, pch=15)# plot diagnostic for q-values
B6 <- eQTLB6[, "adjustedpvalue"]
BXD <- eQTLBXD[, "adjustedpvalue"]
hist(B6, breaks=50, col=rgb(1,0,0,0.5), xlab="adjusted p-values",
ylab="frequency", main=paste("Diagnostic adjusted p-values plot", tissue, condition), las=1)
hist(BXD, breaks=50, col=rgb(0,0,1,0.5), add=TRUE)
# add legend
legend("topright", legend=c("B6","BXD"), col=c(rgb(1,0,0,0.5), rgb(0,0,1,0.5)), pt.cex=2, pch=15)# percentage of significant eQTLs over expressed genes
length(which(eQTLB6[, "adjustedpvalue"]<0.05)) / length(eQTLB6[, "adjustedpvalue"]) * 100## [1] 22.0459
## [1] 22.24566
qB6 <- qvalue::qvalue(eQTLB6$pvalue_secondpermutation)
qBXD <- qvalue::qvalue(eQTLBXD$pvalue_secondpermutation)
summary(qB6)##
## Call:
## qvalue::qvalue(p = eQTLB6$pvalue_secondpermutation)
##
## pi0: 0.6395766
##
## Cumulative number of significant calls:
##
## <1e-04 <0.001 <0.01 <0.025 <0.05 <0.1 <1
## p-value 1361 1790 2540 3049 3566 4348 12819
## q-value 1119 1460 2015 2392 2834 3407 12819
## local FDR 905 1109 1475 1704 1903 2191 12794
##
## Call:
## qvalue::qvalue(p = eQTLBXD$pvalue_secondpermutation)
##
## pi0: 0.6250782
##
## Cumulative number of significant calls:
##
## <1e-04 <0.001 <0.01 <0.025 <0.05 <0.1 <1
## p-value 1364 1785 2538 3053 3547 4349 12764
## q-value 1119 1444 2032 2398 2847 3418 12764
## local FDR 893 1110 1467 1709 1919 2207 12747
# load mm9 cis-eQTL
ceLNSD_mm9 <- read.table("F:/BXD/BACKUP/mm9/Data/IntermediateLayer/cis-eQTL/Liver_nsd/ciseQTL.Liver.NSD.pvalcorrected.txt",sep='')
ceLSD_mm9 <- read.table("F:/BXD/BACKUP/mm9/Data/IntermediateLayer/cis-eQTL/Liver_sd/ciseQTL.Liver.SD.pvalcorrected.txt",sep='')
ceCNSD_mm9 <- read.table("F:/BXD/BACKUP/mm9/Data/IntermediateLayer/cis-eQTL/Cortex_nsd/ciseQTL.Cortex.NSD.pvalcorrected.txt",sep='')
ceCSD_mm9 <- read.table("F:/BXD/BACKUP/mm9/Data/IntermediateLayer/cis-eQTL/Cortex_sd/ciseQTL.Cortex.SD.pvalcorrected.txt",sep='')
# load mm10 cis-eQTL
ceLNSD_mm10 <- read.table("F:/BXD/BACKUP/mm10/Data/IntermediateLayer/cis-eQTL/ciseQTL.Liver.NSD.pvalcorrected.txt",sep='')
ceLSD_mm10 <- read.table("F:/BXD/BACKUP/mm10/Data/IntermediateLayer/cis-eQTL/ciseQTL.Liver.SD.pvalcorrected.txt",sep='')
ceCNSD_mm10 <- read.table("F:/BXD/BACKUP/mm10/Data/IntermediateLayer/cis-eQTL/ciseQTL.Cortex.NSD.pvalcorrected.txt",sep='')
ceCSD_mm10 <- read.table("F:/BXD/BACKUP/mm10/Data/IntermediateLayer/cis-eQTL/ciseQTL.Cortex.SD.pvalcorrected.txt",sep='')
for(data in c("ceCNSD_mm9", "ceCSD_mm9", "ceLNSD_mm9", "ceLSD_mm9", "ceCNSD_mm10", "ceCSD_mm10", "ceLNSD_mm10", "ceLSD_mm10")){
d <- get(data)
print(data)
print(nrow(d))
print(nrow(d[d$adjustedpvalue<0.05,]))
print(nrow(d[d$adjustedpvalue<0.05,]) / nrow(d) * 100)
}## [1] "ceCNSD_mm9"
## [1] 14889
## [1] 4522
## [1] 30.37142
## [1] "ceCSD_mm9"
## [1] 14889
## [1] 4732
## [1] 31.78185
## [1] "ceLNSD_mm9"
## [1] 14103
## [1] 3155
## [1] 22.37113
## [1] "ceLSD_mm9"
## [1] 14103
## [1] 2654
## [1] 18.81869
## [1] "ceCNSD_mm10"
## [1] 15734
## [1] 4192
## [1] 26.64294
## [1] "ceCSD_mm10"
## [1] 15734
## [1] 4542
## [1] 28.86742
## [1] "ceLNSD_mm10"
## [1] 12647
## [1] 3092
## [1] 24.44849
## [1] "ceLSD_mm10"
## [1] 12647
## [1] 2695
## [1] 21.3094
# plot variant-gene expression pair
plotGeneExpressionVariant(ref="B6", gene=rownames(top3genesCNSD), variant=top3genesCNSD$variantB6, tissue="Cortex", condition="NSD")## $B6
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## $BXD
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
# plot variant-gene count pair
plotGeneCountsVariant(ref="B6", gene=rownames(top3genesCNSD), variant=top3genesCNSD$variantB6, tissue="Cortex", condition="NSD")## $B6
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## $BXD
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
# plot variant-gene expression pair
plotGeneExpressionVariant(ref="B6", gene=rownames(top3genesCSD), variant=top3genesCSD$variantB6, tissue="Cortex", condition="SD")## $B6
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## $BXD
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
# plot variant-gene count pair
plotGeneCountsVariant(ref="B6", gene=rownames(top3genesCSD), variant=top3genesCSD$variantB6, tissue="Cortex", condition="SD")## $B6
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## $BXD
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
About 2% of the genes are highly affected. Knowing we affected less than 1% of bp in the genome.
legend: black=B6 allele, tan=D2 allele
Checking for raw gene counts. There are comparable since we used the same reads on B6 or BXD references.
# plot variant-gene expression pair
plotGeneExpressionVariant(ref="B6", gene=rownames(top3genesLNSD), variant=top3genesLNSD$variantB6, tissue="Liver", condition="NSD")## $B6
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## $BXD
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
# plot variant-gene count pair
plotGeneCountsVariant(ref="B6", gene=rownames(top3genesLNSD), variant=top3genesLNSD$variantB6, tissue="Liver", condition="NSD")## $B6
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## $BXD
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
# plot variant-gene expression pair
plotGeneExpressionVariant(ref="B6", gene=rownames(top3genesLSD), variant=top3genesLSD$variantB6, tissue="Liver", condition="SD")## $B6
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## $BXD
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
# plot variant-gene count pair
plotGeneCountsVariant(ref="B6", gene=rownames(top3genesLSD), variant=top3genesLSD$variantB6, tissue="Liver", condition="SD")## $B6
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## $BXD
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
##
## $<NA>
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33
## R version 3.5.3 (2019-03-11)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18363)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=French_Switzerland.1252 LC_CTYPE=French_Switzerland.1252
## [3] LC_MONETARY=French_Switzerland.1252 LC_NUMERIC=C
## [5] LC_TIME=French_Switzerland.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] vioplot_0.3.5 zoo_1.8-7 sm_2.2-5.6 knitr_1.30
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.5 compiler_3.5.3 pillar_1.4.6 plyr_1.8.6
## [5] tools_3.5.3 digest_0.6.25 jsonlite_1.7.1 evaluate_0.14
## [9] lifecycle_0.2.0 tibble_3.0.1 gtable_0.3.0 lattice_0.20-38
## [13] pkgconfig_2.0.3 rlang_0.4.7 yaml_2.2.1 xfun_0.18
## [17] stringr_1.4.0 dplyr_1.0.2 generics_0.0.2 vctrs_0.3.4
## [21] grid_3.5.3 tidyselect_1.1.0 qvalue_2.12.0 glue_1.4.2
## [25] R6_2.4.1 tcltk_3.5.3 rmarkdown_2.4 farver_2.0.3
## [29] ggplot2_3.3.2 purrr_0.3.4 reshape2_1.4.4 magrittr_1.5
## [33] scales_1.1.1 ellipsis_0.3.1 htmltools_0.5.0 splines_3.5.3
## [37] colorspace_1.4-1 labeling_0.3 stringi_1.4.6 munsell_0.5.0
## [41] crayon_1.3.4